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ABSTRACT 

This paper reviews the simulation of earthquake occurx’ence by 
numerical and laboratory mechanical block models. Simple 
linear rheologies are used with elastic forces driving tlie main 
events and 'idscoelastic fbrees being important for aftershock and 
creep occurrence. Friction and its dependence on velocity, stress, 
and displacement also plays a key role in determining how, when, 
and where feult motion occurs. The discussion of the qualitative 
behavior of the simulators focuses on the manner in which energy 
is stored in the system and released by the unstable and stable 
sliding processes. The numerical results emphasize the statis- 
tics of earthquake occurrence and the correlations among source 
parameters. 
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NUMERICAL AND LABORATORY SIMULATION 
OE FAULT MOTION AND 
EARTHQUAICE OCCURRENCE 

L INTRODUCTION 

This paper reviews numerical and experimental simulations of fault motion 
and earthquake occurrence. In the simulations discussed here fault regions are 
represented by discrete blocks which are driven across a friction surface by 



various combinations of elastic and viscous force elements. This representation 
of earthquake sliding and associated fault dynamics lias both strengths and 
weaknesses whiohi can be stated from the outset. Foremost among its strengths 
are its clarity and conceptual simplicity. The modeling of &.ult dynamics in 
terms of sliding blocks is readily accomplished both in the laboratory and on a 
computer. The details of the sliding process appeal to intuition and experience 
and thereby permit a clear understanding of the observed behavior. Many of the 
large scale phenomena associated with earthquakes and fault motion are mini- 
iced in the simulations including the occurrence of foreshocl?s, main shocks , and 
aftershocks in single earthquake sequences, the correlation among source pa- 
rameters such as that between energy drop and the product of rupture length 
and average displacement, and the occurrence of stable, aseismic, sliding epi- 
sodes. In addition, the simulation technique provides a convenient test bed for 
examining new models of earthquake source mechanisms and for exploring 
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tlieir dynamic consequBUces. Tlxe simplicity of the technique also accounts for 
its wealmesses. Aside from the obvious dangers of ovex’simplification of com- 
plicated tectonic and geologic processes, the lumpinoss of tlie block elements 
introduces restrictions of tlie length scale of the phenomena which can be studied 
and introduces a discretization wluch may not be present in nature. 

Willie this review will concentrate on the mechanical block representation 
of earthquake occimrence, note should be made of tlie dislocation theory tech- 
niques for studying some of tlie same problems that will be discussed hex’e. 
Dislocation techniques usually, although not necessarily, assume some Imown 
ruptime propagation characteristics along tlie surface expression of a fault and 
deduce the consequential motions of points elsewhere. Dislocation techniques 
have been discussed by Haskell (1969) among others. 

The organization of tliis paper is as follows. In the remainder of the in- 
tx’oduction we present a simple matliematleal analysis of the sliding occui’ring 
in an eartliquake, the model consistiiig of a single mass block being driven by 
elastic loading. This model serves as a basis for flie more sophisticated mod- 
els intx'oduced latex’. With tlie basic ideas introduced by this model in mind, 
we begin in Section H a systematic discussion of simulation concepts. This 
section discusses linear x’lieological models of the stx’ess-straiii relations in 
x’ocks. Particular emphasis is place on the time dependent behaviox*. Section 
H also reviews various concepts mid experiments oonoerning the fx'ictionai 
x’esistance to sliding, particulax’ly stiok-sllp sliding which appears to be 
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important in earthquake dynamics. In Section n we discuss in a qualitative 
manner the belmvior of various simulators which use the rheological models 
discussed in the previous section. Then in Section IV we discuss simulator 
statistics and correlations among source parameters as deduced from the ex- 
periments and calculations, where possible the results are compared with 
results for naturally occurring events. In the final two sections we discuss 
several topics related to simulator analysis, present our conclusions, and 
comments on the possibilities for future research. 

We begin the discussion of earthquake fault motion by considering the 
model, shown in Figure 1, of a block sitting on a friction surface and subject 
to stress due to the accumulation of elastic strain within a spring. The spring, 
with elastic constant k, is stretched at constant velocity u. In the ab- 
sence of any initial motion in the block (mass=m), the frictional resistance 
prevents sliding until tlie stress rises to the maximum frictional strength, 
f . Once sliding has begun, the frictional resistance drops to the dynamic 
value, f Thus with the initial position of the block at x s 0, motion is initi- 
ated at time tjj when 

kittp = f® (1) 

and tlience 

rax = kCut-x)-f‘^ (2) 
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Tlie solution of eciuntion 2 is 


X = — ^ — : sin- 

k 


^rrt‘-<o' ^ /IT , rr 

■y— -y— +uu-t„)-uyYsm^- »-t„) (S) 


FOX' practical puiposes, once the sHding event is initiated we can ignore terms 
in u since tiiese terms contx'lbute negligible displacements during the lime span 
of an earthquake (typically u = 5 cm/yri tins approximation is matliematically 
equivalent to setting t « tp on tlxe right side of equation 2 wluch then becomes 
mx + lex = f® - f Hence 


X 



(t-tiO 


From tills equation we deduce tile followingj 
a. duration of the event 


(4) 



b. total displacement 


Ax 


2jP-f^) 

k 


(5) 


( 6 ) 


0 , I’eduotion in dx'iving tension, AF, or fx’aotioual driidng tension, AF/Fq 



d. potential energy drop 

AE = P'Ax = 


(S) 


e. peal: velocity 


*nu\ 


f. peak acceleration 


ma\ 


Ilrii 

Vink 


P ~ f 


m 


( 9 ) 


( 10 ) 


Notice tliat tlie reduction in the driving spring tension is independent of the 
elastic constant and the block mass; in fact, the mass affects only the duration 
of the event, tlie velocity and acceleration. Illustrating the relationsliips by 
talcing m a S n 10^** gm, k » 1 x 10^® dyne/cm, f * = 2 x lO’® dyne, f® = 0.8 f% 


(corresponding roughly to a ten kilometer cube at shallow deptlis) 

At ~ 5 sec 
Ax = 40 cm 

AF = 4 X 10^*’ dyne 
AF 


F 


= 0.4 


AK = 3 X I o'* ei'gs 
^nax ~ llcavsec 
= 20cui/sec“ 


This model of blocking sliding under the influence of spring loading intro- 
duces many of the concepts advanced in the more complicated models discussed 
later. Itor now, however, we turn to a discussion of the fundamental stress- 
strain relations for rocks. 
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H. LINEAR BHEOLOGIC MODELS OF STRESS AM) STRAIN IN ROCKS; 
FRICTION 

The Diij'd.eling of feult motion requires as a firndameEtal starting point some 
relationsMp between ihe stress, or, aotiiig on the rocks and the consequential 
strain, e*. In general these relationships can be quite complicated but con- 
siderable understanding can be achieved by restrictiug the analysis to relation- 
ships which are linear, that is those which involve stress, strain, or their time 
derivatives to only the first power. Mathematically this is an extremely power- 
ful assumption for it permits either the total stress or strain due to several ele- 
ments to be determined by summing individual contributions. The relationships 
that we state will be for the case of uniaxial compression, no greater generali- 
zation will be required at present The simplest case, that of pure elasticity 
is, 

a = ke (11) 

A substance obeying this relationship is said to be linear elastic or Eookean. 

A linear viscous substance, or Newtonian viscous, is defined by 

cr = ??6 (12) 

More complicated substances can. be constructed from combinatious of these 
elements, and the very powerful linear operator techniques can be used to de- 
rive the stress-strain relations. We write an element stress-strain relation 

*In accordance with common practice we will often mix the concepts of force and stress. In the present context the 
diffciencc between tire tensor nature of stress and the vector nature of force is not of great signifleonec. 


as 


a - Oe (13) 

where O is a linear operator on e. For the elastic element Oe = k/ edt = 

Icle = ke, while for the viscous element Oe ** ije. In formal manipulations the 
inverse of the integration operator is the differentiation operator; i. e. , 

ID = DI = 1. For elements connected in series as in Figure 2a, the strains, 
a."d hence the strain rates, add directly, i.e., e =. e,. + or 

e = [(kl)“i + u (I4av 

Hence 

a +— ff = ke (14b) 

n 

This is the stress-strain relation for the Maxwell viscous substance formed 
from the series combination of Hookean and Newtonian elements. As shown in 
Figure 3, the strain Increases linearly with time when the stress is held con- 
stant in such a substance. The stress decreases exponentially with time con- 
stant t = tj/Ic when the strain is held constant 

When elements are connected in parallel as in Figure 2b ihe stresses add 
directly. Thus for this so called Kelvin element 

o = 7?e + ke (15) 

In Contrast to the behavior of the Maxwell substance, the Kelvin substance under 
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constant stress obeys the relation 

eg = eoe-‘/^+YO-e-‘/^); =Y ( 16 ) 

Thus in the absence of any stress the initial strain relaxes exponentially while, 
if the substance is initially unstrained, the strain rises toward an asympotic 
value a/k (Figure 3). When the strain is held constant the stress is the same as 
ihat for the linear elastic element. The Kelvin element is the simplest model 
to give transient creep shown, for example, by the time dependent terms in 
Equation 16 (Jaeger and Cook, 1976). It is not adequate however for represent~ 
ing seismic activity because it does not permit instantaneous elastic strain 
(Jaeger and Cook, 1976). 

There are two non-degenerate three element substances. These are shown 
in Figure 2c and 2d. It should be noted that there are alternative equivalent 
representations of these models. For example, Figure 2c can be replaced by 
a spring, kj, in parallel with a series combination of a spring, k^, and a dash- 
pot (viscous element) provided 



The stress-strain relationship for the generalized Kelvin substance is 
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Tims imdex’ conditions of constant stress 


e = ej,e- 




(18) 


and under constant strain 


J ss 


Ooe 


■t/r' 


+ JL 

K K 


d_e“t/r’); / = 


rt 


K + K 


(19) 


Examining equations 18 and 19 we see that this substance shows transient creep 
and transient stress although for t > r, t ' both the stress and strain approach 
steady state values (Figure 3). Instantaneous stress and strain changes are also 
possible. Thus this model ip useful in discussing transient creep and time 
dependent stress readjustment following major sliding events. These features 
will prove to be important in modeling aseismic fault motion and aftershocks. 

For the other non~degenerate three element substance 


. . I' ’'.”1. .. , ’1.'= 

Q g - ^ g ^ 


’lb ^b 


( 20 ) 


Thus for constant stress 


6 = 6o + 


-~+ Aq - — V(I -e"'/’’) 


; T =■ 


% 


( 21 ) 


and for constant strain 


a = CqC**/’’; t = 




k 


( 22 ) 


wliich has tlie same fonn as the Maxwell substance under constant strain. 
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Four element moclola onn also bo used to simixlato features of fault dynamics. 
Figure 2e sliows one form, tlie Burgers substauce, whiob. consists of a series 
combination of Koivin and Maxwell elements. The steess-strabi relation is 


a + — 


« ^ “ or 


ni, \ 


... , . 

0 = k„e + k„ e 


This mode], sliows instantaneous strain and both transient and steady state creep 
but because of its complexity has not been used in a detailed simulation of tlie 
type currently mider disoussion. 

The linear operators we have used here are directly onalngous to tliose ap- 
pearing in odier fields, partioularljf electrical circuit tlieory. The equation 
0 = Oe is anaingous to die electrical equation J ^ ctE where J is the oieetric cur- 
rent density and E is the oieetric field amplitude. Thus stress is analagous to 
current density and stuatn is analagous to the electric field amplitude. The 
operator O is analagous to conductivity or the electrical admittance. Its inverse 
is resistivity or Impedance. Thus die viscous element widi a direct proportion- 
ality between a and e plays die role of a resistor, and die spring tliat of an in- 
duotive coil. It can also be shown dint the inertial effect of die block mass, m, 
is analogous to die capacitive effects hi electrical cixmiits. 

In die proceeding discussion we have developed equations relating stress 
and strain for a munbor of model substances. If die mechanical block dis- 
cussed ill Sootiou I were fx-ee to move widiout encouutering fx^ctional resistance 
then it would undergo a continual slow movoment to provont any long term 
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stress accumulation. This would not, however, be an adequate model of fault 
surfeoes where frictional interactions between the opposite sides of the fault act 
to restrict motion. Thus we need to include in the simulator a description of 



the frictional properties of rocks in the feult zone. We again consider a block 
model with the block resting on a siu’face and subject to a compressive normal 
stress, c. When a shear stress is applied parallel to the contact plane, it is 
found Qiat over a wide range of normal stresses and temperatures and for a wide 
range of rock materials, that no motion occurs until some critical magnitude of 
shear stress, F*, is reached. An experimental relationship between Fs and a 
can be deduced viz, 

F® = (24) 

where p® is the coeffloient of static friction. Expei'iments have found that to a 
good approximation p® is only a wealc function of tlie area of contact and the nor- 
mal stress (Jaeger and Cook, 1976). Once sliding is initiated, it is again ob- 
served that tile frictional resistance is proportional to normal stress but, as a 
general rule, the proportionality constant, p** , Icnown as the dynamic friction 
coefficient is less than the static coefficient. The value of the dynamic friction 
coefficient may vary with tb.e speed of sliding. 

When careful measurements are made of the relationship between static 
fi'iction and normal stress it is found that the F® versus or cimve does not pass 
through the origin, but that to a better approximation F® = Fq + where F^ 
corresponds to an intrinsic shear strength. 
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‘ ^ tion. In. general one miglit expect that the static frictiou coeffiLoient win increase 


p with the duration of contact following some earlier sliding, corresponding to a 

1 time dependent healing of the ruptured ragiotu Dieterich (1972a) has observed 

I 1 . ■ 

this effect and concludes 

I M®(t) = p®(0) [1 +alag(At)J; At> 1 sec. (25) 

where At is the duration of contact following sliding, 

i ' 

I A different approach to the frictional resistance in block sliding has been 

I discussed by Byerlee (1970), Nur and Shultz (1973), and Nur (1977) in light of 

f 

• ! T ■ ■ ■ 

earthquake dynamics. In this approach it is assumed that asperities along the 
fault siu:faoa lock together to form the frictional resistance. When the shear 
force becomes strong enough to cause brittle failure of these asperities motion 
can occur. As the sHdiiig proceeds, continual locking and brealdng of the asperi- 

i : ■ 

> ties cause fluctuations in the frictional resistance. Thus, the frictioiml resist- 

I 

ance varies in some irregular manner with displacement, m this model the 
value of the frictional resistance at the initiation of the unstable sliding event is 
analogous to the static friction and the average value of the fluctuating friction 
during the earthquake sliding corresponds to the dynamic friction. An interest- 
ing feature of ihis model is that both rapid, unstable sliding resembling that oc- 
curing in earthqualces and slow, stable sliding resembling diat in creep events 
J can occur. The criterion for unstable sliding is that the drivii^ stress decease 

less rapidly with displacement than the frictional resistance. 

In the earthquake simulations, various models for the frictional resistance 

i will be employed. In the simplest eases, static and dynamic friction forces, f® 

'll 







and f'\ lire used to explain the gross features of the main sliding events. Other 
models are introduood for modeling aftershooks, creep, and the details of the 
carthciuake instability mechanisms. || 

in. MODEL SIMUXATIOMS - GENERAL BESULTS 

The prototype for the models that are discussed in this paper is shown in 
Figure 4. It consists of a number of mosses oomieoted together by springs and 
resting on a friction surlaco. The masses are driven by spriiigs comiected to 
a moving plate. In some eax’ly laboratory simulatious only the first mass of die 
olioin was connected to a driving source. In many of die computational models 
die springs have been replaced by more complicated rheological elements such 
as the viscoelastic elements disoussod in die previous section. Details concern- 
ing die various laboratory and computational models tiiat are discussed in this 
review can be found in the Appendix, 

The simplest models of fiiult motion during eardiquidces employ elastic 
springs as coupling elements and frictioiial surfiices over wliich the blocks 
slide. Laboratory models of dlls typo have been studied by Burridge and 
Knopoff (1967), King and Kuopoff (19G8b), and King (1975). Computational 
models include those of Dieterioh (1972b, 1978), Otsuka (1972a), Cohen (1977a, 

1977b), and Ruudle and Jackson (1977). In examining the sequence of events 
occurring in their simulator, Burridge mid Kuopoff (19G7) found that small events 
occur largely at random while large events involving major changes in elastic 
potential energy occur nearly periodically. Between major events die poten- 
tial energy of the system increased nearlj' linearly widi time when die hiter- 
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mass springs Iiad equal elastic constants. Furfiieimore, the major events all 
occurred at approximately i3ie same level of potential energy and involved 

comparable energy drops. Ilius the period between these large events was 
simply the energy drop divided by the rate of energy input. 

King and Knopoff (1968b) also found in their study that the strain energy 
density before a shock showed little variation with the magnitude of the event. 
An analagous result has been reported by Bath and Duda (1964) for natural 
events. King (197S) observed a long term linear trend in average displacement 
with time; this behavior is, of course, necessary since in the long term the 
displacement rate must keep pace with velocity of the moving splate. More 
significantly, he noted that large shocks occurred at times when the displace- 
ments of most masses had fallen below their long term averages. On the other 
hand King, unlike Burridge and Knopoff (1967), reported that large shocks were 
r Dt periodic in time. 

Using a computational model in which friction was allowed to vary from 
block to block, Dieterich (1972b) found that the frequency with which blocks 
move is related to their static frictional strength with low friction blocks 
moving more frequently than those with hi^ friction. * Dieterich (1973) also 
perform^dd two and three dimensional simulations and derived scaling laws that 
relate source parameters obtained with one set of input variables to source 
parameters with a new set of input parameters. These important scaling laws 
are summarized in Table I. 

^Onc coiitmst between the laboratory and computational models should be noted. In the laboratory mo deb stick-slip 
sliding is observed but there Is no certainty about the matlicmatical dcsciiption of tlic friction law. In the computational 
model the friction law is clearly defined mathematically} but it is not certain how well the mathematical description 
represents the natural process. 
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The importance of seismic gaps, that is regions of low seismic activity in 
an otherwise highly seismic region, has long been studied for eaarfihqualce pre- 


diction purposes. Otsuka (1972) found that large simulation eairthquakes occurred 
in regions of anomalously low seismic activity. 

Studies have also been performed to determine how the qualitative chamc- 
teristics of feult movements are influenced by the degree of heterogeneity in the 
friction and elastic parameters. Cohen (1977b) has found, for ejjample, that 
with relatively homogeneous (ault parameters successive events tend to propa- 
gate along the length of the fatilt in a type of epicenter migration. He also found 
that recurrence patterns in the locations and magnitudes of seismic activity 
could be observed. Both effects were less common with heterogeneous faults, 
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but seismic gaps were more likely to occur. 

Bundle and Jackson (1977) searched for evidence of an increase in simula- 
tor seismicity preceeding large events such as su^ested by studies of Oallforaia 
sefemicity by Wesson and Ellsworth (1973) and Wyss and Lee (1973). Such be- 
havior was not present in their elastic model leading them to speculate that the 
seismicity increases may be due to (among other causes) time dependent fault 
parameters. 

Several workers, including Bundle and Jackson (1977) have pointed out that 
as the ratio of elastic modulus to stress increases, there are frequent small 
displacements and the system behavior becomes ductile. Conversely as the 
ratio decreases larger infrequent events occur. 
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Ill the models that have been discussed up to now the only significant time 
dependent processes have been the slow elastic loading of the driving springs 
and the sudden elastic motions diming simulation earthquakes. As a consequence 
the models are not capable of ejqilaining phenomena occuring on intermediate 
time scales. I3ie most important of these intermediate scale phenomena are 
earthquake sequences involving either foreshocks, the main shock, and after- 
shocks or nearly equal magnitude swarms. Other important time dependent 
phenomena include continuous, episodic, premonitory, and post seismic creep. 
To study some of these eSects, various computational models of viscoelasticity 
and time and stress dependent friction have been employed. lu the earliest of 
these models Burridge and Knopoff (1967) divided their one dimensional fault into 
three regions, two regions with predominately elastic properties and low viscos- 
ity separated by a region of high viscosity. The effect of the viscous region was 
to introduce a time delay between the time of sliding in one elastic region and 
the time at which the resultant stress change was felt in the other elastic re- 
gion. This time delay was responsible for the occurrence of aftershocks both 
in ihe secondary and in the primary elastic regions. In the high viscosity re- 
gion, unstable sliding could not occur. The time dependent characteristics of 
the energy release in the viscoelastic model were significantly different from 
the elastic case. In a sequence of 69 aftershocks following a main shock on 
the primaiy jEault, tiie authors reported; 

1. The cumulative energy released by aftershocks increased as (1 - e“*^’). 
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2. The total stress energy released in the shock sequence was divided as 
follows: main shock - 12%, aftershoclcs on primary segment - 44%, 
aftershocks on secondary segment - 13%» Additionally a potential 
energy release equal to 43% of the energy x’eleaSed by simulation 
shocks was dissipated by viscous flow, 

3. Despite the greater energy release in the primary segment, only 
27% of the shocics, including the main shock, occurred along this 
segment. 

4. The decay in the frequency of aftershocks occurring with time was 
fairly consistent with either tiie (1 + jSt)'^ dependence suggested by 
Omoil (1894) or an exponential decay. 

It may be significant to note that in this sequence more potential energy is re- 
leased by aftershocks than by the main shock. We will comment more on this 


point later. 


Dieterich (1972b) has proposed that the interplay of viscoelastic stress x’e- 


oovery following an earthqualce coupled with a reduced frictionEil strength of tlie 
fault following a sudden sliding may be responsible for aftershocks. Consider 
die generalized Kelvin substance discussed in section H. . It can be shown that 
following a sudden stress relaxing sliding event, a time dependent Increase in 
stress can occur with the stress rise asymptotically approaching k^ / + kj,) 

of the stress drop during the sliding. Furthexmore, as Dietei'ich's (1972a) 
experimental work suggests, the fidctional resistance against further sliding 




is less liian the resistance pilor to sliding. Presumably this is related to the 
finite time required for frictional healing of the rupture (equation 25). if 
the viscoelastic stress recovery is more rapid than the time required for fault 
healing an aftershock may be generated. The aftershocks occur in the same 
region of the fault that ruptured during the original event and the duration of 
the aftershock sequence increases with the size of the main shock. Further- 
more a correlation has been found between the dispacements and magnitudes 
in the main shock and those of the largest aftershock's (Cohen, 1977b). 
Although the frequency of aftershocls decreases with time, it does not follow 
Omoi’i's law (Dieterich, 1974). 

hi an alternative model of aftershocks, Bundle and Jackson (1977) focused 
attention on stress induced crack nucleation. They assumed that the friction 



strength, F®, against sliding decreases at a rate proportional to the amount 
the stress exceeds some value of friction, F®“ , i.e., 


dF®(t) 

dt 


1 


k(t)-F®“] 

r 


(26) 


where F®“ is a constant and o(t) is the stress acting on the block. In this 
model, if the dynamic friction is greater than F®“ , the stress drop during the 
first rupture is small and the stress, cr(t), remains large causing a decrease in 
frictional strength. Since friction is then reduced while stress remains high 
aftershoclcs may occur. The aftershocks are not confined to the region of the 
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adjacent to that which ruptured during the primary event (this border region 
having been stressed by the elastic elements coupling it to the displaced region). 
Since the stress, a (t), generally changes little on the time scale of t, t is also 
the time scale for ihe aftershock occurrence. Bundle and Jackson^s results 
show too few aftershocks to examine die dependence of aftershock frequency on 
time following the main shock, however, it may be necessary to introduce some 
time variation in <T(t), other dian the tectonic loading rate; to produce a deca3dng 
frequency of aftershoclcs resembling that present in nature. 

IV. STATISTICS, SOUKCE PABAMETEES, AND COERELATIONS 


Up to this point our discussion of simulator results has been primaidly 
descriptive. Quantitative results generally fell into two categories - one, the 
statistics of earthquake occurrence and two, the correlations among source 
parameters. Considering first the frequency of simulator earthquake occur- 
rence as a function of the energy released in the events Burridge and Khopoff 
(1967) obtained the experimental results shown in Figure 5. Except for the 
lowest energies the data fit a straight line wiih 


log f = A — B log E 


(27) 


where f is the frequency of events occurring with energy greater than or equal 
to E. This relationship appears analagous to the relationship between the fi'e- 
quency of naturally occurring events and the seismic energy log f = A'- b' log E^. 
Bearing in mind the relationship between seismic energy and earthqualce 
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m.E.i.’iiitude, viz, log Ej = a + bM it is tempting to define a magnitude for the 
simulator events by 

^^SIM = log E ( 28 ) 

This analogy suggests a Imear relationship between log E and log E^ and while 
this is a useful viewpoint it is ondy an approximation whose validity needs to be 
confirmed. A reason for the flattening of the simulator frequency versus mag- 
nitude curve at low magnitudes is that the discrete nature of the model imposes 
some lowest magnitude event below which no events will occur. 

For naturaUy occurring events B' = 0. 4 while Burridge and Knopoff (1967) 
found B = 1 in their mechanical simulator. These authors have conjectured 
that the difference between their simulator results and actual earthquake data 
may be a reflection of the one-dimensional nature of their model. However 
for a very similar mechanical model, King and Knopoff (1968) found B = 0. 71 
and for a later laboratory model King (1975) found B = 0. 5. Bundle and Jackson 
(1971) concluded ihat the results from their computational simulator could be 
fit to two lines wilh B = 0. 1 - 0. 5 for low magnitude events and B = 1-5 for high 
maguitude events. An alternative analysis by Otsulca (1972b) revealed domward 
curvature in the logarithm of frequency versus magnitude curve, the degree of 
downward curvature decreasing with the probabiliiy that the motion of one blocks 
stimulates an adjacent block into motion. A similar curvature was found in the 
computational simulator by Cohen (1977a). Although it appears that there is 


I * 

i : 


some qualitative agreement between the deviations from a linear beliavior for 
some actual earthqualce data and the numerical results, the quantitative agree- 
ment is at best only feir. 

The variation in the frequency of the time intervals between events can be 
examined to determine whether the distribution of interevent periods is 
Poissionian as would be expected if the events occur randomly in time. Both 
Burridge and Knopoff {1967) and Bundle and Jackson have found some non- 

Poissonian components to the distribution, presumably due to the interactions 
between adjacent events. The motion of one block alters the stretch or com- 
pression of the spring connecting it to the neighboring block thus altering the 
stress and the time at which the stress overcomes the fricilonal resistance. 
However, ICnopoff, Mitchel, and Jackson (1972) and Bundle and Jaclcson (1977) 


f i 


Imve shown tliat the occurrence rate for the simulation events are consistent 
with a stoclmstic model operating on the stored elastic energy. They have used 


in their analyses tlie Kolmogorov backward equation 


Ji(E)P(E) + a = f*''"^"'‘'T(E,x)X(x)P(x)dx 

rIF •'n 


- E 


where now E is the elastic potential energy of the system, X(E) is the probability 


that an event will occur in time interval dt if the enei^ is E, a is the late at 
which enei'gy accumulates between events, T(E,x) dE is the probability that if 


an event occurs at initial energy the final energy will be in the range E to 
E + dE and P(E)dE is the probability of finding the system energy between 
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E and E + dE. Figure 6 shows the degree with which the Kolmogorov equation 


fits the results of Eundle and Jackson (1977). 

King and Khopoff (1968b) used their mechanical simulator to deduce corre- 
lations among the rupture parameters for model shoclcs. Shock energy is 
shown in Figure 7 as a function of the length of tne rupture, £, or the number of 
displaced masses, L. The data shov/s a concave up behavior on the log-log 
plot of Figure 7; however, if the data with L = 8 are excluded a linear correla- 
tion with log E = 1. 3 + 1. 6 log L can be used. Since L = 8 corresponds to 
movement within the entire spring-mass system and since the free boundary 
conditions of the end masses are not representative of natural conditions, the 
data witli N = 8 must be viewed with caution. Dieterich (1974) has assembled 
a number of ms^nitude versus rupture length relations used by various 
workers. In general they have the form m = g + h log £ with h between one 
and two. 

Somewhat better linear correlations are discovered when log E is plotted 
against log LD, log LD^, and log Ld where D and d are the pealc and average 
displacements of the blocks. The correlation between E and Ld is shown in 
Figure 8. The relationship E « £d will be established below while the relation 
E a i’ (or E a £ d^) has been deduced from calculations based on an electro- 
static analog (Knopoff, 1958). In a study of earthquake data King and Knopoff 
(1968a) deduced 


logXD^ = 2.24M-4.99 


(30a) 
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wlillo in liieir simulator work (King and Knopoff, 19GSb) tlioy found 

logLD^ = 1.55 log 1.89 (SOb) 

Assimiing log ICjj = 11. 8 *t* 1. 5M and EwtjEj wlioro ij is tlio soiamic efficiency, 
til© simulator results givo 

logLD^ = 2.5M+ 1.64 -h 1.55 log 4. (80o) 


Thus the slopes of tlie curves of log LD- versus M are comparable for tlie 
natural and simulated events provided the assumed rolationahip between, 
seismic energy and simulator euiergy is valid. 

Some of tlie results presented In die px'ocoodlng paragraph can be explained 
by relatively simple tlieorotloal arguments. For simulations in whioh the elastic 
energy stored in die dxdviiig and connecting springs is dissipated by friction 
sliding, tlio energy drop during an event is 



From this we can deduce that E« £D“; here 1’ is the dlmonslon of the 
block. This result agrees witli flie aforemontionod result by linopoff (1968). 
In the second case we assume the friction varies little from block to block so 
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it can be removed from under the summation sign and 


E = f^'EAx, = f^Ld (33) 

a result agreeing witli the correlation suggested by ITigure 8, 

The preoeeding analysis can also be used to explain the simulation observed 
variation in average displacement versus magnitude. From equation 33 we 
deduce 


d = — , a Jlii « constant (34) 

f‘*L L 

Alternatively for single block motion 
^ IT 1o“sim/2 

d a-A/ - oc one block (35) 

’ VF 

The simulation data shown in figure 9 are consistent with tliese equations. 

The change in die elastic force in the spring-mass systems due to the un- 
stable sliding in an event is simply 

AF = X; kf Axj (36) 

i 

or if the didving spring moduli are all approximately equal 

AF = k“Ld (37) 


Thus in the case where both equations 34 and 37 are valid AF «= AE = 10^’®*'''’ while 
for single block motion AF « 10 ^'sim^? One observed variation in the stress drop 
in the driving springs versus log E is shown in Figure 10, In genei'al the data 
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for L < 8 suggest that while the stress drop increases with m^nitude it does 



so more slowly than indicated by the preoeeding equations. 

Different simulations have indicated somewhat different variation in 
displacement with rupture length. Results obtained by both King and 
Knopoff (1068b) and King (1975) revealed both increasing average and peak 
displacement with increasing rupture length. By contrast Cohen (1977b) 



discovered tliat the average displacement increased with rupture length 
for small numbers of displaced blocks but for longer ruptures found no 
further increase in average displacement. Wallace (1974) has assembled 
data relating anipture length and maximum displacement to earthquake magni- 
tude for earthquakes in Western North America. His data indicates that 
maximum displacement increases about linearly with rt5)ture length although 
there is considerable data scatter. These results agree better with those of 
King (1975) and King and Knopoff (1968b), than those of Cohen (1977b). 

V. OBSERVATION AND DISCUSSION OF RELATED TOPICS 

It is interesting to consider for a moment the mathematiail origin of the 
instability that characterizes the sliding process in various simulations. In the 
simplest models, the instability is due to the sudden transition from the static 
to the lower dynamic friction when the driving force equals the static friction, 
Burridge and Knopoff (1967) used a more complicated friction law and deduced 
that the condition for instability on a single block is that the net friction force 
decrease with increasing velocity at the beginning of unstable sliding. 
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(A similar friction law wMcIi is mathematically smoother than that of Burridge 
and Knopoff but which reduces to their law in the limits of small and large veloc- 
ities has been used by "Whithead and Gans (1974) ). Thus the simple model is con- 
tained within the more general treatment of Burridge and Knopoff. Notice that 
the instability condition focuses on the variation in stress for a change in veloc- 
ity. A somewhat different view advanced by Nur and coworkers and alluded to 
in Sectioa n, focuses on how the stress changes with displacement. If F(x) is 
the net diriving force and f(x) a position dependent friction then the condition for 
unstable sliding in a region of decreasing friction is -f'(x) > -F'{x). For the purely 
elastic case in which F(x) = fc(ut - x), this condition is -f'(x) > k. This insta- 
bility condition implies, in the elastic case, that friction decrease more 
rapidly with displacement than the reduction in the elastic spring driving tension. 
The most detailed quantitive caiculatiott in Nur's so called stif&iess model has 
been perfoinned with a one block viscoelastic substance (Cohen, 1978). The 
most interesting features that emerge are 

1. creep occurs premonitory to unstable sliding episodes. Transient 
creep may also occur after unstable sliding. 

2. stable sliding occurs on both epidsoic and long-term time time scales. 
Episodic sliding events may be stable or unstable depending on the 
interplay between the viscoelastic and friction parameters. 

Despite these apparent correlations between the viscoelastic stiffiiess model 
and naturally occurring events, it must be emphasized that there is 




i 
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considerable uncertainty about tile basic validity of the stiffiiess theory. Mod- 
el calculations can only deduce the consequences of the assiunptions and it 
can happen that alternailve assumptions have the same consequences. 

ICnopoff (1972) advanced a tlieory of aftershock occurrence which contains 
some of the elements of other aftershock models and in a sense anticipated 
some of the later developments. The two key elements in KnopofPs model are 
the occurrence of large stresses on die edges of a nature zone of a shock and 
the delay in rupture of a highly stressed region by prerupture tertiary creep, 
in application KhopofE notes that because of the elastic properties the boundary 
between a region that slid during aai earthquake and the neighboring unmoved 
region is one of Mgh stress concentration. Knopoff, like Handle and Jackson 
(1977), postulated that fliis Mglily stressed boundary region may not fail in- 
stantaneously but may fail witli a delay time dependent on the degree of over- 
stressing. This delayed response gives rise to aftershocks which cluster 
around die edges of the main shock rupture region. As already noted failure, 
m this model, occurs with preshock tertiary creep. 

A common featui'e of the viscoelastic aftershock theories that we liave re- 
viewed in tiiis paper is the occurrence of either creep or stress recovery fol- 
lowing the primary shock. The fundamental study of the relationship between 



aftershocks and creep using the linear circuit theory that we presented in Sec- 
tion II is due to Benioff (1951). In applying his theory to a number of earthquake 


sequences Benioff concluded that wliile the seismic wave energy of the primary 
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shock exceeds that of the entire aftershock sequence, the elastic strain energy 
dissipated in the creep-aftershock sequence may be comparable to that of the 
primary shock. In the earlier discussion of the numerical model of Burridge 
and Knopoff (1967) we pointed out that more potential energy was released in 
the aftershock and creeping following the main shock than was released in 
the main shock. These two results strongly suggest that aftershocks and post- 
seimic creep may be responsible for a much larger fraction of the total elastic 
energy dissipation than would be suggested by their seismic wave magnitudes. 

As we mentioned earlier the spring and block models that have been re- 
viewed in this paper are just a small subset of all models of fault motion and 
earthquake occurrence. As one alternative approach rock mechanics experiments 
involving the sliding of small rock samples have provided a wealth of data on 
the rupture and sliding processes. In the last several years Brune (1973) has 
used stressed foam rihber to observe both stick slip sliding and creep. Along 
computational lines, dislocation theory has been used to determine ground 
displacements, velocities, and accelerations due to fault motions. The models 
discussed here and dislocation theory models become somewhat similar in con- 
cept when the present models are e3q>anded to two and three dimensions and the 
elements made small, and when the dislocation models are formulated with a 
deseriptiou of the rapture process. 
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VI. CONCLUSIONS 

It has been somewhat over a decade since the first massive block simulator 
of earthquake occurrence was developed, hi that time well over a dozen studies 
have been made by almost as many researchers. Major contributions have been 
made in understanding the propagation of rupture along a feult, in developing 
correlations among source parameters, and in explaining the interplay between 
foreshocks, main shocks, and aftershocks and creep. It appears to the author 



that at least two areas of study have potential for contributing significant new 
information. The first of these is the study of two and three dimensional visco- 
elastic systems, possibly with multiple faults. It seems likely that the interac- 
tion of elements not directly on the foult surface will produce a number of ef- 
fects such as a propagating creep wave not observed in most cue dimensional 
simulations. The second area is the study of the physics connecting die near 
field block sliding with the far field seismic wave motion. A synthesis of tech- 
niques such as a combination of a sliding block description of rupture with a 
dislocation theory description of die far field effects might prove fruitful. By 
contrast it seems less lilcely diat fundamentally new insights will be achieved 
by utilizing more complicated rheological models (althou^ simple, but novel 
models may be illuminating). The more complicated the rheological model, 
die greater the number of adjustable parameters that impact the calculations. 
Thus entirely different models may produce similar effects if the parameters 
are appropriately adjusted, hi this case, the validity of one model as opposed 
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to another cannot be established. Finally it seems increasingly evident that a 


detailed understanding of earthquake occurrence will require improvements in 



the present Icnowledge of friction laws, rupture criteria, and the physics of 
highly stressed rock. 
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APPENDIX 


The various laboratory and numerical simulations discussed in this paper 
are summarized below. 


Burridge and Knopoff (1967) 
a. laboratory model 


One dimensional - 8 masses connecling by springs and driven by 
stretching of spring connected to first mass; m = 142 gms; spring 
constant = 2 x lo® dynes/strain; spring lengths: 3 cm or 1.5 - 12 cm; 
loading rate: 2 cm/sec. 
b. computational model 

One dimensional - 10 masses connected by springs and driven by 
springs coupling to a moving plate. 5 blocks represent a primary 


fe-ult, and 3 blocks represent an auxiliary fault. The primary and 
auxiliary fault segments are separated by tv/o blocks in a highly viscous 


region, hi a normalized ssrstem of un.ts m = 1, connecting spring 
constants = 100, driving spring constants 1, intermass spacing = 1. 
The friction law is 


f(x) = < - 


E X 


1+A(x-H) 


1 - A(x + H) 


IxKH 


x>H 


x<-H 




with A = 10, E = 1 (all blocks); H = 10"® (1-3, 6-10), 10^ (4, 5) 

B = 5 (1-3), 10^ (4, 5), 15(6), 10 (7-10), when the numbers in paren- 
theses are the block numbers for which the quoted values are appli- 
cable. The driving rate in this system of units is 10“® . 

2. Kings and Knopoff (1968b) 

Laboratory model similar to that used by Burridge and Knopoff (1967) ; 8 
masses: m = 130 gms; spring constants = 10^ dynes/cm; spring un- 
stretched length = 3 cm. 

3. Otsuka (1972) 

Two dimensional computational model; blocks suspended to fixed overhead 
support by leaf springs and in contact with moving floor. 2000 blocks in a 
100 X 20 grid. 

4. Dieterich (1972b) 

One dimensional computational model; 50 bloclcs interconnected and driven 
by coupling to moving plate. Four models are employed in which the 
coupling elements and the friction vary as follows : 

a. elastic springs, time independent static friction 

b. elastic springs, time dependent friction 

c. viscoelastic elements, time independent friction 

d. viscoelastic elements, time dependent friction. 

The viscoelastic coupling element is the generalized Kelvin element. The 
time dependent friction is represented as f® (t) = f®(tg) [1 + A log t] for 
t > 1 sec. 
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The numerical values used in the model calculation are: friction stress 
- 100-300 bars, = 1. 5 - 5. 0, t = 10^ - 10“^ sec, = 0. 8f® for 

time independent friction, t= 0. 98 f® (tg) for time dependent friction, 
A - 0.02, interelement spacing = 1 Ion, driving plate velocity = 5 cm/yr 
Dieterich (1973) 

Two and three dimensional computational models in which the boundary 
conditions are specified displacements of bloclcs far removed from the 
fault Elastic springs are assumed, length to widtii ratio of the I’egion 
varies from 2-4. 

King (1975) 

One dimensional laboratory model in which masses are driven by coupling 
to an overhead rotating flat circular plate; 8 masses, m = 110-113 gm; 
driving spring constant « 2 x 10** dyne/cm; connecting spring constant « 

1. 6 X 10** dyne/cm. Connecting springs under an average tension of 
1. 3 X 10* dynes, driving plate velociiy = 0. 03 cm/sec. 

Cohen (1977b) 

One dimensional computation model similar to Dieterich's (1972) models a. 
and d. Spring constants = 10 - 10 * ® dynes/cm; 1. S^kjAi <t < 

lO'* sec; 0 <f® < 3 X 1020 dynes; 0 <fVf®< . 99 
Cohen (1977a) 

One dimensional computational models using elastic spi’ings- Spring con- 
stants S X 10*® dynes/cm - 1. 5 x 10*”* dynes/cm; f® = 1 x 10^® dynes - 
3 X 10^® dynes; f** = 0. 8 f*. 
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9. Bundle and Jadcson (1977) 

One dimensional computational model employing 20 massless elements 

di’iveu by moving plate both an elastic and unelastic model are invested. 

a. elastic model: driving spring constants = 100-10, 000 kbars/cm; con- 
necting spring constant 0-S200 Idb/cm; static friction 0. 3 Icb - 9 Id). 

b. anelastic model: The frictional strength is time and stress dependent 
with 

dps 1 

— [o(t) 

at r 

where a(t) = time dependent stress on the block, and F®“ and r are 
constant. 

In the choice of numerical values for the spring constants in tlie various 
simulations, little attention has been devoted to the relationsliip between the 
driving and connecting spring constants. However, Yaniashlta (1976) iias ana- 
Ized the relationslups between the spring constants and the elasticity constants 
IX and X where (i is rigidity and X the usual Lame paranxeter. For a one 
dimensional model he finds 



whei'e Vp/Vj. is tlie ratio of pi’iniary to secondary seismic wave velocities, and 
Az and Ay are the block dimensions along and perpendicular to ilie fault. The 
parameter Az is tlie dimension of the block in tiie third direction. 
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Table L 


Scale factors 

lengib 

a = l/l' 

stress 

cr 

II 

elastic constants 

c = ju/ju' = 

density 

d 

Scaling relations 

force 

F = a^bF' 

spring constants 

k = ack' 

displacement 

Ax = abc"^ Ax' 

mass 

m = a^dm' 

velociiy 

Ax = b(cd)"^Ax 

acceleration 

Ax = b(ad)-^ Ax 

energy 

E = a^b-c-iE' 

seismic moment 

M = a^bM' 


Scaling laws for source parameters. Source 
parameters are Imown in pilmed system and 


desired in unprimed system (from Dieterieh, 19731 






Burgers substance 


Figure 2. Multiple element rheologic models: a. Maxwell substance, b. Kelvin 
substance, c. three element generalized Kelvin substance, d. alterna 
tive three element substance, e. four element Burgers substance 
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Figure 6: Plot of the values of the Kolmogorov backward equation versus nondimensional energy E, where E 
has been broken up into 10 energy bands of equal width. The top curve is the left side of the 
equation and the bottom curve is the residual, defined as the left side minus the ri^t side. (Bundle 
and Jackson, 1977) 
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Figure 8: Shock enex’gy versus the 
dlsplacemeut. Solid dots 
and Knopoff, 1968) 
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